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Abstract: The application of finite-difference 
time-domain methods to the analysis of structures 
having curved boundaries is a long-standing 
problem. Traditionally, either the staircasing 
approximation which is simple but inaccurate, or 
the generalised nonorthogonal grids which are 


accurate but complex have been used. In this 
contribution a novel approach is presented which 
combines the efficiency of the Cartesian mesh 
with the accuracy of the conformal grid. Results 
are presented for some example structures which 
are in good agreement with other methods. 


1. ‘Introduction 


Finite-difference time-domain (FDTD) methods origi- 
nally put forward by Yee [1] have proved to be very 
efficient numerical algorithms in computational electro- 
magnetics. However, the traditional FDTD algorithm 
is based on a Cartesian co-ordinate system, and it is 
difficult to generate meshes exactly for electromagnetic 
structures with curved boundaries. It is usual to utilise 
a staircase approximation in the FDTD method for 
curved structures and an accurate solution can only be 
obtained by using very fine grids and consequently, a 
very small time step. In addition, the staircase approxi- 
mation of the physical boundary of a resonator will 
often result in a failure to detect all the resonant modes 
and in the prediction of Q factors which are lower than 
in reality. Apart from the staircase approximation, 
there are two main modifications to the FDTD method 
which have been put forward to analyse the elec- 
tromagnetic structures with curved boundaries: the 
contour path (CPFDTD) algorithm [2-5], and the non- 
orthogonal FDTD method [6-9]. The CPFDTD algo- 
rithm is based on the integral form of Ampere’s and 
Faraday’s laws. The update equations must be modi- 
fied in those distorted cells which are near the curved 
boundaries. This is readily achieved for curved perfect 
conductors because the tangential components of the 
electric field on the boundaries are zero. Even so, 
although good results have been obtained by Jurgens, 
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[2, 3], the generation of the mesh is nontrivial and, in 
addition, it has been shown that the basic scheme must 
be modified to ensure stability [4, 5]. For dielectric 
objects, where the tangential £ field is nonzero, addi- 
tional equations are needed to calculate the tangential 
components of electric field from the equation of 
Ampere’s law. This makes the algorithm much more 
complex and unwieldy and, with the exception of a 
simple 2D example in [2], no published results are 
known. to the authors. 

The generalised nonorthogonal FDTD algorithm was 
put forward by Holland [6] and refined by others [7, 8]. 
The underlying mathematics dealing with the geometry 
and its application to electromagnetic fields can be 
found in [10]. Many numerical results have been 
obtained successfully [7, 8, 11]. However, compared 
with the conventional FDTD scheme, two additional 
equations are needed in each iteration step to realise 
the transform between the contravariant and covariant 
components of electric and magnetic fields. Moreover, 
extra computer memory is needed to store the metric 
tensor of both E and H nodes, which are essential 
parameters in nonorthogonal co-ordinate systems. For 
these reasons, the method is limited to relatively small 
structures. Although the nonorthogonal FDTD scheme 
is a generalised method, it is difficult to make it com- 
patible with existing FDTD software which is based on 
the Cartesian co-ordinate system. 

To take advantage of the nonorthogonal FDTD 
algorithm which is powerful for the computation of 
structures with curved boundaries, we have developed a 
combined method, which uses the theory of the non- 
orthogonal FDTD scheme within an underlying Carte- 
sian co-ordinate system. In the combined method most 
of the grid is in the Cartesian co-ordinate system with 
only those cells near the curved boundaries being’ 
treated as nonorthogonal cells. The relative advantages 
and disadvantages of the new method compared with 
existing techniques are given in Table 1. 


Table 1: Advantages (+) and disadvantages (-) of. 
methods 


CPFDTD Nonorthogona! This research 
Computer 

+ - + 
resources 
Accuracy 7 + + 


General material 


: + + 
boundaries 


Ease of mesh 
generation 
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The fundamental concepts and some simulation work 
using the combined method on closed curved metal 
structures have been presented in [12]. In this contribu- 
tion the method is developed to allow the electromag- 
netic analysis of dielectrically loaded waveguides with 
arbitrary cross-sections, which include dielectric rod 
and slabs inside rectangular and circular waveguides. 
The simulation results will be compared with those 
from analytical method and the FDTD methods using 
the staircase approximation. 


2 Description of new method 


The new FDTD algorithm which combines the advan- 
tages of the nonorthogonal FDTD program with con- 
ventional Cartesian FDTD technique has been used to 
compute the cut-off frequencies of circular, rectangular 
and elliptical hollow waveguides. In the combined algo- 
rithm, the majority of the mesh is Cartesian which 
allows the use of the conventional FDTD and only 
near the oblique surfaces are non-Cartesian cells used. 
This leads to a more accurate formulation than 
CPFDTD but much less computational resources than 
a generalised nonorthogonal method. 

A good understanding of the algorithm can be 
started by considering the metric tensor, which is one 
of the important components in the analysis of vectors 
in a nonorthogonal co-ordinate system. It can be 
obtained from g; = 4, - Ay, where A,, A; are the bases 
of a covariant vector. For an orthogonal co-ordinate 
system only the diagonal components will be nonzero 
(as shown in eqn. 1). Moreover, the contravariant and 
covariant components will be collinear in an orthogo- 
nal co-ordinate system. Extra memory is not needed to 
store the metric tensors and contravariant or covariant 
components. Only the values of those parts of the mesh 
which are close to the material boundaries need to be 
stored. 


9:3 = Ave A; 


= |aas] - dij 


_fl isj 
sua {h 153 “ 
6, is the angle between two generalised axes, which it is 
easy to obtain when these axes are straight lines, other- 


wise, the angle can be obtained from their tangential 
lines. 


= |a;a,;| - cos 6;; 


if A;, 4; are orthogonal 


2.1 Generation of the meshes 

Let an arbitrary curve pass through a standard FDTD 
cell so that the original cell is bisected. The two basic 
types of cells near the curved boundary are defined as 
either ‘flag’ or ‘triangle’. In ‘flag’ cells the original cell 
is split by the material boundary into two cells which 
are still quadrilateral. In this case the neighbouring cell 
is extended as shown in Fig. 1 and the EZ, or EZ, node is 
replaced by an E, node on the material boundary. The 
‘second type of cell, called a ‘triangle’, can not be 
treated directly as an FDTD cell, the mesh needs to be 
reconstructed so that all the cells are quadrilateral. To 
do this, an additional point is defined on the material 
boundary. such. as C’ in Fig. 2 and the edges of the 
intercepted cell and its neighbours are modified as also 
shown in Fig. 2. The new nonorthogonal FDTD 
meshes on the curved boundary for those triangular 
and quadrilateral meshes are then obtained. It is not 
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difficult to establish the functions of the two types of 
cells. The ‘triangle’ cell is always placed where the cur- 
vature of a curved boundary changes abruptly. To 
some extent, it looks like the triangular element near 
the curved boundary in the finite-element method. The 
difference is that in FDTD the ‘triangle’ is actually a 
generalised quadrilateral cell, and the vector of the 
electromagnetic field is denoted in a nonorthogonal co- 
ordinate system. The ‘flag’ cells always appear on the 
curved boundary while the curvature changes slowly. 
In general, these two cells are enough to denote a 
curved or an oblique boundary in FDTD computations 
in the underlying Cartesian co-ordinate system. 


Ez 


Fig.1 ‘Flag’ cell extending to neighbouring cells 


Fig.2 


(i, k) (i+1,k) 


Fig.3 Oblique structure with a corner meshed by new FDTD method 


2.2 Treatment of corners in angular objects 

If we consider oblique structures with corners, such as 
the rotated rectangular waveguide shown as Fig. 3, the 
cell containing the corner will be treated as a ‘triangle’, 
but in this case the extra node point is placed precisely 
on the corner. By this means, the exact position of the 
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corners has been taken into consideration, leading to 
improved accuracy. The example of an oblique rectan- 
gular disk modelled by CPFDTD is shown in Fig. 4. It 
shows that the small quadrilateral area labelled A has 
not been taken into consideration in the CPFDTD 
mesh. This will introduce an inaccuracy if the mesh is 
coarse, which is avoided in the method described here. 


Fig.4 Oblique structure with a corner meshed by CPFDTD 


2.3 Choosing the position of the nodes in 
modified cells 

Once the new co-ordinates including those ‘triangle’ 
and ‘flag’ cells are determined, it is easy to set the E 
and H nodes in each cell. Following Yee’s scheme, the 
middle points of the cell boundaries are chosen on the 
x axis and the z axis, respectively, as E,. and the cen- 
troid of the cell as H nodes, but x, y, z cannot now be 
regarded as Cartesian co-ordinates. This arrangement 
of field nodes allows existing FDTD algorithms which 
have been developed for a Cartesian co-ordinate system 
to be immediately incorporated into new FDTD pro- 
gram. These include absorbing boundary conditions 
and near-far field transformations. A further advantage 
of the new FDTD scheme is that it reduces to conven- 
tional FDTD codes when those structures to be studied 
are simple rectangles. Those components of electric 
field (for TE modes) and magnetic field (for TM 
modes) on the interface of material are moved along 
the boundary, which is used to approximate the curva- 
ture of structures in practical computer programming. 
Having determined the node position is of the electric 
and magnetic fields, the components of the metric 
tensor g, in nonorthogonal cells can be calculated and 
stored (shown as eqn. 1). 


3 FDTD iteration formulae in new meshes 


In the previous Section, the meshing of a curved 
boundary on a Cartesian co-ordinate system was 
defined. This approach ensures that the new mesh con- 
forms both to the material boundaries and to the sur- 
rounding Cartesian mesh. Thus the conventional 
FDTD approach can be used on the Cartesian cells 
and nonorthogonal method on the nonorthogonal cells, 
and the values are passed directly to the common 
nodes without the need for interpolation. 

As an example, consider the structure shown in 
Fig. 3. Cell @, k) is a Cartesian node and the tradi- 
tional FDTD iteration formula shown in eqn. 2 can be 
used to obtain the £ and H components. During itera- 
tion, the covariant and contravariant components need 
not be introduced, and therefore the memory require- 
ments are the same as for the standard: FDTD 
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EE* (4,8) 
= ER (i,k) — (Hg G,&) — HFG, &- 1) 
Fett oh) 
= EEG.) + LHS i — HYG 1,8) 
Hy"? (i,k) 
= Hy 7k) + EN +1,k) — EP (,k)| 

_ —_[E"(é,k +1) — E*(,b)] (2) 


ar 
where 6x and éz are variable space steps given by 


én = risa = nF 
62 = laj41 — | 
The contravariant and covariant components of the £ 


and H fields on the nonorthogonal cells can be defined 
as 


=H.-A, 
Hi=#H. 4 


EK; = E . As 

Ei=- fF. Aj 

2,9 =2@,y Or Zz 
where the covariant vector bases of the vector A(A, 
and A,) are defined along the x and z axes, respec- 
tively, and their dual bases, known as the contravariant 
bases (A* and A’), are defined to be orthogonal to the 


original bases A, and A,, which are shown in Fig. 5 
(2D cases only). 


a Az (zaxis) 


Fig.5 Definition of contravariant and covariant vector bases 


Cell @ + 1, k + 1) is a nonorthogonal node, so the 
nonorthogonal FDTD method is used to find the val- 
ues of the E and H components. The equations for this 
cell are shown as 

E*(@+1,k+1)""! 


eine 


y(i+1,k) 


evs 


Ee(¢+1,k+1)""" 
= E*(i+1,k+1)" 


‘at (j,k+1)+ Ayit+i,k+1)]"* 


= 
Guu 
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—H,(i+1,k+1)]""2 


HYG +1,k+1)"2 


= Hi+1,k+1)th 42h. [Sw 
Cir g 


| AG+2k+1) Eit+l,k+)) 
VGez(t+2,k +1) Vgee(t+1,k +1) 
E,(it1,k+2)  E,(itl,k +1) ~ 
) 


(3) 


In the above equations, the contravariant components 
in cell @, k) need not be transformed into covariant 
components because they are collinear in Cartesian 
cells. In nonorthogonal cells, transformation is neces- 
sary, using eqn. 4. 


E,(t+1,k+1) =Gi2E°+1,k +1) 


Gaz 1 1 
bee ae ae ee es ee 
+ ral (i it kel 


1 
+E? (i+1-5,e+1- 


1 
+E? (i+14 fkei4 


: 

2 

© 

2 

: 1 1 

+4h% (i+1-5.¢+143) 


Hy(t¢+1,k+1) = H¥%44+1,k+1) 


Jax 
Gz = ‘Gxz 4 
Vas (4) 


where x, y, z denote generalised co-ordinates in a non- 
orthogonal co-ordinate system, they will reduce to Car- 
tesian co-ordinates in those Cartesian cells. 

For the conventional FDTD algorithm, the bound- 
ary condition of the material interfaces has automati- 
cally been taken into consideration in eqns. 3 and 4. In 
this paper all outer boundaries of waveguides are 
assumed to be perfect metal, so one can let ¢ > © in 
the formulae. For dielectric materials, one simply 
inputs the value of permittivity, and allows ¢ = 1 if it is 
in a free space. It should be emphasised here that all 
the above equations can be rewritten so as to include 
materials of finite conductivity. In this case the meth- 
ods used in conventional FDTD can be directly applied 
to the new algorithm. 


4 Numerical results 


Several numerical examples of useful dielectric-loaded 
waveguides are presented in this Section. The structures 
considered here are: first, dielectric-slab-loaded 
waveguides with homogeneous and inhomogeneous fill- 
ings; secondly, dielectric-rod-loaded square waveguides; 
and finally, circular waveguides loaded with dielectric 
cylinders. The numerical results obtained using the new 
FDTD algorithm will be compared with data from 
other methods such as the variation-iteration method 
[13], analytical method (14], [15] and FDTD method by 
staircase approximation. 


4.1 Dielectric-slab-loaded waveguides with 
homogeneous or inhomogeneous fillings 

Fig. 6 shows the mesh of the slab-dielectric-loaded 
waveguide, where all the cells which are near to the 
material boundary have been made conformal to the 
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material. The position of the slab in the rectangular 
shield with the dimensions b/a = 0.5, a = 0.3m can be 
described as x; = (a — 1)/2, x» = (a + 1/2, where a and 
6 are the width and height of the rectangular 
waveguide, respectively, and ¢ is the thickness of the 
dielectric slab. As presented in the previous Section, 
this combined method shows the advantage of dealing 
with electromagnetic structures with corners which do 
not coincide with the Cartesian mesh boundary. This 
avoids the constraints involved in choosing the grid 
sizes in the standard or graded grid FDTD algorithm 
which result from ensuring that the mesh coincides 
with the material boundaries. Fig. 7 shows the results 
for the cut-off frequencies of the LSE;,) mode as func- 
tions of ¢ and ¢, of the homogeneous dielectric slab. 


03 
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0.1 


Fig.6 Dielectric-slab-loaded waveguide meshed in Cartesian grids 
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Fig.7 Cut-off frequencies of LSEjq mode in dielectric-slab-loaded 
waveguides 


They show a good agreement with the results by Yu 
and Chu [13]. In this case, the dielectric-loaded 
waveguide is meshed by 14 x 14 cells with the computa- 
tion area of 0.38m x 0.38m. Table 2 shows the cut-off 
frequencies of some higher order modes existing in the 
homogeneous dielectric-slab-loaded waveguide with 
é, = 2. Results by FDTD algorithm with staircase 
approximation are also listed in Table 2. The approxi- 
mate algorithm fails to detect all of the resonant modes 
and shows less accuracy than the combined FDTD. 
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Table 2: Results of cut-off frequencies (GHz) of LSE,,,, modes in a 


homogeneous dielectric centrally loaded waveguide with 


blaz=tla= 1/2, x, = (a- 0/2, x) = (a+ 0/2, and e{x) =2 


Combined Combined Staircase Staircase Staircase From 
Modes FDTD FDTD FDTD FDTD FDTD 113] 

(14 x 14) (28 x 28) (14 x 14) (28 x 28) (50 x 50) 
LSE19 0.3718 0.3718 _— — —_ 0.3701 
LSE59 0.7984 0.8080 0.7348 0.7855 0.7940 0.8058 
LSE,, 0.8151 0.8240 _— — — 0.8224 
LSE; 1.1846 1.1883 1.0929 — — 1.1261 
LSEx39 1.2346 1.2608 4.1144 — — 1.2739 
LSE, 1.4634 1.5039 1.3562 1.3898 1.4981 1.4970 
LSE3, 1.6183 1.5397 1.4939 —_ _— 1.5244 
LSE,o 1.7018 1.6636 1.5835 1.6615 1.6674 1.6978 
LSE. —‘1.8066 1.7542 1.6771 — — 1.7331 
LSE,, 1.8948 1.9091 1.7890 1.8283 —_ 1.9083 


It is interesting to note that in some cases the use of a 
finer staircased mesh leads to more modes being 
missed. The corresponding results in the offset dielec- 
tric-loaded waveguide are listed in Table 3. Both of 
them show a good agreement with values given by [13]. 


Table 3: Results of cut-off frequencies (GHz) of LSE,,,, 
modes in a homogeneous dielectric offset loaded 
waveguide with b/a = t/a = 1/2, x, =f, xX) = a, and €{x) =2 


From combined 


Modes EDTD (28 x 28) From [13] 
LSEqo. 0.3980 0.4029 
LSEo) —-0.8413 0.8489 
LSE. 0.8699 0.8666 
LSE», 1.2012 1.2151 
LSE9 1.2251 1.2195 
LSE34 1.4396 1.4837 
LSE12 1.4849 1.5291 
LSEgq 1.6755 1.6750 
LSE29 1.8185 1.8315 
LSE,, 1.8972 1.8691 
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Fig.8 Plot of cut-off frequencies of LSE,g mode in an inhomogenous 
dielectric-slab-loaded waveguide 


Calculations were also carried out for the case of an 
inhomogeneous dielectric filling with ¢, = 4(€,max — 1) 
(x — x1)(X2 ~— x)/(x2 — x1). The cut-off frequencies of the 
waveguide of LSE), mode with various slab thickness 
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and different dielectric ¢, ma, values were calculated. In 
the simulation, it is difficult to obtain more precise 
results for those thinner slabs in the waveguide, as 
many more sampling points are needed to describe the 
variations of the inhomogeneous dielectric material. 
The results for the cut-off frequencies showing an aver- 
age error of approximately 1.5% for a 14 x 14 mesh 
and 0.5% for a 28 x 28 mesh on waveguides with inho- 
mogeneous dielectric fillings are presented in Figs. 8 
and 9. 
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4.2 Dielectric-rod-loaded square waveguides 

A square waveguide (side length A) with a centred die- 
lectric rod (radius ag) will be investigated in this sec-' 
tion. This exemplifies the meshing of a dielectric-loaded 
electromagnetic structure with curved boundaries 
within a Cartesian mesh. The mesh is shown in Fig. 10, 
which indicates that the computational region of 0.38m' 
x 0.38m is meshed by 14 x 14 cells. The modes in this. 
structure are very complicated and their properties 
have been studied in [14]. The simulation results using 
the combined FDTD method are displayed in Fig. 11 

where they compared with values from [14] and also 
show good agreement with them. All expected higher- 
order modes are detected. The displayed results show 
the cut-off frequencies of HE modes in the guides plot- 
ted against rod radius and confirm the validity of the 
combined FDTD algorithm. 
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Fig.10 Dielectric-rod-loaded square waveguide structure meshed in Car- 
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Fig.11 Cut-off frequencies of different modes in the dielectric-rod- 

loaded square waveguide 
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4.3 Circular waveguides with dielectric 
cylinders 

Fig. 12 shows the structure of a uniform circular guide 
consisting of two concentric dielectric cylinders. The 
radius of the inner one is a with a relative permittivity 
€,, and the radius of the outer one is } with a relative 
permittivity e. The dominant modes in it are angularly 


symmetric Eo, or Ho, modes [15]. To determine the cut- 
off frequencies it is possible to solve the transcendental 
equations which contain radial functions. This is a 
complicated process compared with numerical meth- 
ods. Table 4 gives the cut-off frequencies of the domi- 
nant TM mode in the waveguide calculated using the 
combined algorithm and the staircase approximation of 
the FDTD method using two different cell sizes. This 
shows that with the staircase approximation, the solu- 
tions only converge when a fine mesh is used whereas, 
with the combined method accurate results are 
obtained using only 12 x 12 mesh. 


Fig. 12 Circular guide with dielectric cylinders meshed in Cartesian 
gri 


5 Conclusion 


A new FDTD scheme has been introduced, which is 
based on using a nonorthogonal FDTD algorithm so 
as to deal with an arbitrary electromagnetic structure 
on an underlying Cartesian grid. The validity has been 
confirmed by the simulation of dielectrically-loaded 
waveguides with arbitrary cross-sections. The new com- 
bined method has the benefit of reduced computer 
resources and easily generated meshes when compared 
with schemes based on generalised nonorthogonal co- 
ordinate systems and of improved accuracy when com- 
pared to a staircase approximation. The numerical sim- 
ulation shows that the new method is very efficient and 
the results are in excellent agreement with those from 
other methods. The application and extension of the 
method to three dimensional structures is being studied 
and will be reported in a future contribution. 


6 Acknowledgment 


The authors wish to thank the Committee of Vice- 
Chancellors and Principals (UK) for the provision of 


Table 4: Cut-off frequencies (GHz) of the lowest E mode in 
circular guide with dielectric cylinders (b/a = 2, €, = 1) 


& Combined Combined 
FDTD (12 x 12) FDTD (14 x 14) 

2.54 0.5822 0.5912 

10.2 0.3195 0.3231 
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53.0 0.1430 0.1447 

96.0 0.1073 0.1070 
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FDTD staircase 
approximation 


FDTD staircase 
approximation 


(12 x 12) (72 x 72) 
0.5510 0.5938 
0.3027 0.3357 
0.2050 0.2247 
0.1360 0.1491 
0.1012 0.1069 
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